Jonathan, Environmental Data Science @UIUC. Lab 10 - Emission point analysis
4/5/2024
Instructions
Go through this R file in detail and understand the meaning of each script line. Answer the four questions (Q1 to Q4). Then organize the entire R file into a R Markdown file. Your R Markdown file should include each of the steps as below, including your answers to the four questions.
# set working directory
setwd("/Users/jonathanbernard/Desktop/UIUC/Spring 2024/GGIS 224/GGIS224/Labs/Lab10")
# load libraries
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(sf)## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(tmap)## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')library(leaflet)# COUNTY BOUNDARY POLYGON DATA
# 2018 Illinois County Boundaries
# Source: Census, https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
counties.original <- st_read("Data/County Boundaries/cb_2018_us_county_500k.shp")## Reading layer `cb_2018_us_county_500k' from data source
## `/Users/jonathanbernard/Desktop/UIUC/Spring 2024/GGIS 224/GGIS224/Labs/Lab10/Data/County Boundaries/cb_2018_us_county_500k.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3233 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.1489 ymin: -14.5487 xmax: 179.7785 ymax: 71.36516
## Geodetic CRS: NAD83head(counties.original) # KEY = GEOID is 5-digit county ID ## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -89.18251 ymin: 36.91554 xmax: -83.45285 ymax: 38.52538
## Geodetic CRS: NAD83
## STATEFP COUNTYFP COUNTYNS AFFGEOID GEOID NAME LSAD ALAND
## 1 21 007 00516850 0500000US21007 21007 Ballard 06 639387454
## 2 21 017 00516855 0500000US21017 21017 Bourbon 06 750439351
## 3 21 031 00516862 0500000US21031 21031 Butler 06 1103571974
## 4 21 065 00516879 0500000US21065 21065 Estill 06 655509930
## 5 21 069 00516881 0500000US21069 21069 Fleming 06 902727151
## 6 21 093 00516893 0500000US21093 21093 Hardin 06 1614569777
## AWATER geometry
## 1 69473325 MULTIPOLYGON (((-89.18137 3...
## 2 4829777 MULTIPOLYGON (((-84.44266 3...
## 3 13943044 MULTIPOLYGON (((-86.94486 3...
## 4 6516335 MULTIPOLYGON (((-84.12662 3...
## 5 7182793 MULTIPOLYGON (((-83.98428 3...
## 6 17463238 MULTIPOLYGON (((-86.27756 3...# Subset to just IL Data
counties.IL <- counties.original |> filter(STATEFP == 17)
plot(counties.IL)
# ILLINOIS STATE BOUNDARY
il <- st_union(counties.IL)
plot(il)
# COUNTY HEALTH ATTRIBUTE DATA
# 2023 Length of Life Estimates
# Source: County Health Rankings, https://www.countyhealthrankings.org/explore-health-rankings/illinois/data-and-resources
# Technical Documentation: https://www.countyhealthrankings.org/sites/default/files/media/document/2023%20CHRR%20Technical%20Document.pdf
# Length of Life RANK = lower scores indicate BEST health, higher scores indicate WORSE health
lengthLife <- read.csv("Data/IL_LengthLife_CHR.csv")
head(lengthLife) #FIPS = County Key## FIPS State County LengthLife LLRank
## 1 17001 Illinois Adams -0.10 46
## 2 17003 Illinois Alexander 1.50 102
## 3 17005 Illinois Bond -0.13 40
## 4 17007 Illinois Boone -0.40 18
## 5 17009 Illinois Brown -0.80 4
## 6 17011 Illinois Bureau -0.16 38# Merge Health Data to master county file
counties <- merge(counties.IL, lengthLife, by.x="GEOID", by.y="FIPS") #Essential a table join.
head(counties) #Didn't work; inspect...CO_FIPS a bit strange!## Simple feature collection with 6 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -91.51308 ymin: 36.9703 xmax: -88.70541 ymax: 42.49505
## Geodetic CRS: NAD83
## GEOID STATEFP COUNTYFP COUNTYNS AFFGEOID NAME LSAD ALAND
## 1 17001 17 001 00424202 0500000US17001 Adams 06 2214804824
## 2 17003 17 003 00424203 0500000US17003 Alexander 06 609996948
## 3 17005 17 005 00424204 0500000US17005 Bond 06 985073312
## 4 17007 17 007 00424205 0500000US17007 Boone 06 727171389
## 5 17009 17 009 00424206 0500000US17009 Brown 06 791704469
## 6 17011 17 011 00424207 0500000US17011 Bureau 06 2250916543
## AWATER State County LengthLife LLRank geometry
## 1 41767689 Illinois Adams -0.10 46 MULTIPOLYGON (((-91.51297 4...
## 2 44237171 Illinois Alexander 1.50 102 MULTIPOLYGON (((-89.51839 3...
## 3 6462629 Illinois Bond -0.13 40 MULTIPOLYGON (((-89.63926 3...
## 4 3361806 Illinois Boone -0.40 18 MULTIPOLYGON (((-88.94098 4...
## 5 4139310 Illinois Brown -0.80 4 MULTIPOLYGON (((-90.91703 3...
## 6 11606402 Illinois Bureau -0.16 38 MULTIPOLYGON (((-89.86235 4...# Quick map for confirmation
# Length of Life RANK = lower scores indicate BETTER health, higher scores indicate WORSE health
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties) + tm_polygons("LLRank", palette = "BuPu") +
tm_layout(frame = FALSE, legend.outside = TRUE)
# ADDING EMISSIONS POINT DATA
# 2020 National Emissions Inventory, Facility Data in IL
# Source: EPA, https://www.epa.gov/air-emissions-inventories/2020-national-emissions-inventory-nei-data
emissions.original <- read.csv("Data/NEI2020_FacilityIL.csv")
head(emissions.original)## pointID SITE..NAME EIS.Facility.ID STATE
## 1 1 ADM Animal Nutrition 3344211 Illinois
## 2 2 Prince Agri Products Inc 3345811 Illinois
## 3 3 Prince Minerals Inc 7314711 Illinois
## 4 4 Quincy Regional Airport - Baldwin Field 3346911 Illinois
## 5 5 Bunge North America Inc 3347011 Illinois
## 6 6 Cairo Dry Kilns Inc 3362711 Illinois
## State.County Pollutant Pollutant.Type Emissions..Tons.
## 1 IL - Adams Lead Compounds CAP/HAP 0.00001
## 2 IL - Adams Lead Compounds CAP/HAP 0.00006
## 3 IL - Adams Lead Compounds CAP/HAP 0.00256
## 4 IL - Adams Lead Compounds CAP/HAP 0.03723
## 5 IL - Alexander Lead Compounds CAP/HAP 0.00032
## 6 IL - Alexander Lead Compounds CAP/HAP 0.00002
## Facility.Type Street.Address
## 1 Unspecified 1000 N 30th St
## 2 Unspecified 4618 Gardner Expy
## 3 Unspecified 401 N Prince Plz
## 4 Airport 1645 Highway 104
## 5 Food Products Processing Plant 203 34th St
## 6 Unspecified 2 Miles N Of Cairo on Hwy 51
## NAICS EPA.Region FIPS Latitude Longitude
## 1 Other Animal Food Manufacturing 5 17001 39.94534 -91.36624
## 2 Other Animal Food Manufacturing 5 17001 39.87616 -91.39827
## 3 Synthetic Dye and Pigment Manufacturing 5 17001 39.89511 -91.40916
## 4 Other Airport Operations 5 17001 39.94093 -91.19425
## 5 Soybean and Other Oilseed Processing 5 17003 37.01605 -89.17777
## 6 Home Centers 5 17003 37.04587 -89.18468#Convert CSV to Spatial Data
emissions <- st_as_sf(emissions.original,
coords = c("Longitude", "Latitude"),
crs = 4326)
# Check for type of pollutant in the dataset
unique(emissions$Pollutant)## [1] "Lead Compounds"# Using dots, adjust style, palette, and number of bins (n)
# By now, you are expected to be able to self-interpret what following statement does.
tm_shape(il) + tm_borders(lwd = 2) + #Draw state boundary
tm_shape(counties) + tm_borders(lwd = 0.5) + #overlay county boundaries
tm_shape(emissions) + #overlay point coordinates
tm_dots("Emissions..Tons.", palette = "BrBG", #then apply graduated color based on attribute "Emissions..Tons."
n = 10, style = "quantile",
title="Lead Emissions (tons)") +
tm_layout(frame = FALSE, legend.outside = TRUE)
# Using bubbles, adjust size, color, style, and number of bins (n)
# Again, you are expected to be able to self-interpret what following statement does.
tm_shape(il) + tm_borders(lwd = 2) + #Draw state boundary
tm_shape(counties) + tm_borders(lwd = 0.5) + #overlay county boundaries
tm_shape(emissions) + #overlay point coordinates
tm_bubbles("Emissions..Tons.", col = "tomato1") + #then apply graduated symbols based on attribute "Emissions..Tons."
tm_layout(frame = FALSE, legend.outside = TRUE)
# Add the health data for exploration of associations with lead emissions.
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties) + tm_polygons("LLRank", palette = "Greys") +
tm_shape(emissions) +
tm_bubbles(size = "Emissions..Tons.", col = "Emissions..Tons.",
palette = "PuRd", style = "quantile") +
tm_layout(frame = FALSE, legend.outside = TRUE) 
tmap_mode("view")## tmap mode set to interactive viewing# Be sure to "Zoom" to explore map interactively
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties) + tm_polygons("LLRank", alpha = 0.5, palette = "Greys") +
tm_shape(emissions) +
tm_bubbles(size = "Emissions..Tons.", col = "Emissions..Tons.",
palette = "PuRd", style = "quantile") +
tm_basemap("Esri.WorldTopoMap")## Legend for symbol sizes not available in view mode.# Search for alternate basemaps --
# Use Leaflet Provider Demo at https://leaflet-extras.github.io/leaflet-providers/preview/
# & try out different Provider Names (in the tm_basemap argument)
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties) + tm_polygons("LLRank", alpha = 0.5, palette = "Greys") +
tm_shape(emissions) +
tm_bubbles(size = "Emissions..Tons.", col = "Emissions..Tons.",
palette = "PuRd", style = "quantile") +
tm_basemap("OpenStreetMap.HOT") ## Legend for symbol sizes not available in view mode.# Switch back to plot mode
tmap_mode("plot")## tmap mode set to plottingst_crs(emissions)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]st_crs(counties)## Coordinate Reference System:
## User input: NAD83
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]]emissions.3435 <- st_transform(emissions, 3435) #"NAD83 / Illinois East (ftUS)"
counties.3435 <- st_transform(counties, 3435)
st_crs(emissions.3435)## Coordinate Reference System:
## User input: EPSG:3435
## wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",36.6666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-88.3333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.999975,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",984250,
## LENGTHUNIT["US survey foot",0.304800609601219],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["US survey foot",0.304800609601219],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["US survey foot",0.304800609601219]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["US survey foot",0.304800609601219]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["United States (USA) - Illinois - counties of Boone; Champaign; Clark; Clay; Coles; Cook; Crawford; Cumberland; De Kalb; De Witt; Douglas; Du Page; Edgar; Edwards; Effingham; Fayette; Ford; Franklin; Gallatin; Grundy; Hamilton; Hardin; Iroquois; Jasper; Jefferson; Johnson; Kane; Kankakee; Kendall; La Salle; Lake; Lawrence; Livingston; Macon; Marion; Massac; McHenry; McLean; Moultrie; Piatt; Pope; Richland; Saline; Shelby; Vermilion; Wabash; Wayne; White; Will; Williamson."],
## BBOX[37.06,-89.28,42.5,-87.02]],
## ID["EPSG",3435]]st_crs(counties.3435) #102 counties## Coordinate Reference System:
## User input: EPSG:3435
## wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",36.6666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-88.3333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.999975,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",984250,
## LENGTHUNIT["US survey foot",0.304800609601219],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["US survey foot",0.304800609601219],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["US survey foot",0.304800609601219]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["US survey foot",0.304800609601219]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["United States (USA) - Illinois - counties of Boone; Champaign; Clark; Clay; Coles; Cook; Crawford; Cumberland; De Kalb; De Witt; Douglas; Du Page; Edgar; Edwards; Effingham; Fayette; Ford; Franklin; Gallatin; Grundy; Hamilton; Hardin; Iroquois; Jasper; Jefferson; Johnson; Kane; Kankakee; Kendall; La Salle; Lake; Lawrence; Livingston; Macon; Marion; Massac; McHenry; McLean; Moultrie; Piatt; Pope; Richland; Saline; Shelby; Vermilion; Wabash; Wayne; White; Will; Williamson."],
## BBOX[37.06,-89.28,42.5,-87.02]],
## ID["EPSG",3435]]nrow(emissions.3435) #850 emission facility points## [1] 850nrow(counties.3435) #102 counties## [1] 102# Count all of the facilities that intersect a county.
# we then use lengths() to find out how many items are present in a vector (i.e. county).
# Then, we store the facility counts for each of the counties into a new attribute called TotEmisFac.
counties.3435$TotEmisFac <- lengths(st_intersects(counties.3435, emissions.3435))
# Map the facility count per county as a choropleth map
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties.3435) + tm_polygons(col = "TotEmisFac", n=10,
style = "jenks", title = "Total Facilities") +
tm_layout(frame = FALSE, legend.outside = TRUE)
# Calculate County Area
counties.3435$county_area <- st_area(counties.3435) Q1:
What is the unit of the calculated county area? Answer: the unit is square meters because the st_area() returns area values in square meters by default even though hte CRS used has units in US survey feet, the st_area() function still returns the area in square meters.
Q2:
Make a choropleth map showing the density of emission facilities (i.e., the count of emission points per square mile) by county.
counties.3435$density <- counties.3435$TotEmisFac / (as.numeric(counties.3435$county_area) / 27878400)
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties.3435) + tm_polygons(col = "density", n=10,
style = "jenks", title = "Facility Density (per sq mile)") +
tm_layout(frame = FALSE, legend.outside = TRUE)
# Generate a point-in-polygon (PIP) dataset, joining county data to point data
pip <- st_join(emissions.3435, counties.3435, join = st_intersects)
head(pip)## Simple feature collection with 6 features and 29 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 121147.9 ymin: 128297.8 xmax: 737702 ymax: 1208455
## Projected CRS: NAD83 / Illinois East (ftUS)
## pointID SITE..NAME EIS.Facility.ID STATE
## 1 1 ADM Animal Nutrition 3344211 Illinois
## 2 2 Prince Agri Products Inc 3345811 Illinois
## 3 3 Prince Minerals Inc 7314711 Illinois
## 4 4 Quincy Regional Airport - Baldwin Field 3346911 Illinois
## 5 5 Bunge North America Inc 3347011 Illinois
## 6 6 Cairo Dry Kilns Inc 3362711 Illinois
## State.County Pollutant Pollutant.Type Emissions..Tons.
## 1 IL - Adams Lead Compounds CAP/HAP 0.00001
## 2 IL - Adams Lead Compounds CAP/HAP 0.00006
## 3 IL - Adams Lead Compounds CAP/HAP 0.00256
## 4 IL - Adams Lead Compounds CAP/HAP 0.03723
## 5 IL - Alexander Lead Compounds CAP/HAP 0.00032
## 6 IL - Alexander Lead Compounds CAP/HAP 0.00002
## Facility.Type Street.Address
## 1 Unspecified 1000 N 30th St
## 2 Unspecified 4618 Gardner Expy
## 3 Unspecified 401 N Prince Plz
## 4 Airport 1645 Highway 104
## 5 Food Products Processing Plant 203 34th St
## 6 Unspecified 2 Miles N Of Cairo on Hwy 51
## NAICS EPA.Region FIPS GEOID STATEFP
## 1 Other Animal Food Manufacturing 5 17001 17001 17
## 2 Other Animal Food Manufacturing 5 17001 17001 17
## 3 Synthetic Dye and Pigment Manufacturing 5 17001 17001 17
## 4 Other Airport Operations 5 17001 17001 17
## 5 Soybean and Other Oilseed Processing 5 17003 17003 17
## 6 Home Centers 5 17003 17003 17
## COUNTYFP COUNTYNS AFFGEOID NAME LSAD ALAND AWATER State
## 1 001 00424202 0500000US17001 Adams 06 2214804824 41767689 Illinois
## 2 001 00424202 0500000US17001 Adams 06 2214804824 41767689 Illinois
## 3 001 00424202 0500000US17001 Adams 06 2214804824 41767689 Illinois
## 4 001 00424202 0500000US17001 Adams 06 2214804824 41767689 Illinois
## 5 003 00424203 0500000US17003 Alexander 06 609996948 44237171 Illinois
## 6 003 00424203 0500000US17003 Alexander 06 609996948 44237171 Illinois
## County LengthLife LLRank TotEmisFac county_area
## 1 Adams -0.1 46 4 24325392002 [US_survey_foot^2]
## 2 Adams -0.1 46 4 24325392002 [US_survey_foot^2]
## 3 Adams -0.1 46 4 24325392002 [US_survey_foot^2]
## 4 Adams -0.1 46 4 24325392002 [US_survey_foot^2]
## 5 Alexander 1.5 102 3 7044304829 [US_survey_foot^2]
## 6 Alexander 1.5 102 3 7044304829 [US_survey_foot^2]
## density geometry
## 1 0.004584247 POINT (133815.6 1208455)
## 2 0.004584247 POINT (123966.1 1183559)
## 3 0.004584247 POINT (121147.9 1190566)
## 4 0.004584247 POINT (181996.7 1205255)
## 5 0.011872740 POINT (737702 128297.8)
## 6 0.011872740 POINT (735782.5 139171.5)# Using the PIP file, average all of emissions per county.
temp <- aggregate(pip$Emissions..Tons., by = list(pip$GEOID), mean)
head(temp)## Group.1 x
## 1 17001 9.965000e-03
## 2 17003 6.516667e-03
## 3 17005 1.269750e-02
## 4 17007 2.994800e-02
## 5 17009 3.390000e-03
## 6 17011 2.666667e-05# Rename the data in temp file
names(temp) <- c("GEOID", "EmissionAve")
# Join back to master county dataset
counties.3435 <- merge(counties.3435, temp, by = "GEOID", all.x = TRUE)
head(counties.3435)## Simple feature collection with 6 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95633.61 ymin: 111619.1 xmax: 883631.1 ymax: 2123580
## Projected CRS: NAD83 / Illinois East (ftUS)
## GEOID STATEFP COUNTYFP COUNTYNS AFFGEOID NAME LSAD ALAND
## 1 17001 17 001 00424202 0500000US17001 Adams 06 2214804824
## 2 17003 17 003 00424203 0500000US17003 Alexander 06 609996948
## 3 17005 17 005 00424204 0500000US17005 Bond 06 985073312
## 4 17007 17 007 00424205 0500000US17007 Boone 06 727171389
## 5 17009 17 009 00424206 0500000US17009 Brown 06 791704469
## 6 17011 17 011 00424207 0500000US17011 Bureau 06 2250916543
## AWATER State County LengthLife LLRank TotEmisFac
## 1 41767689 Illinois Adams -0.10 46 4
## 2 44237171 Illinois Alexander 1.50 102 3
## 3 6462629 Illinois Bond -0.13 40 4
## 4 3361806 Illinois Boone -0.40 18 5
## 5 4139310 Illinois Brown -0.80 4 2
## 6 11606402 Illinois Bureau -0.16 38 3
## county_area density EmissionAve
## 1 24325392002 [US_survey_foot^2] 0.004584247 9.965000e-03
## 2 7044304829 [US_survey_foot^2] 0.011872740 6.516667e-03
## 3 10673939103 [US_survey_foot^2] 0.010447277 1.269750e-02
## 4 7864025133 [US_survey_foot^2] 0.017725274 2.994800e-02
## 5 8576601423 [US_survey_foot^2] 0.006501037 3.390000e-03
## 6 24360918399 [US_survey_foot^2] 0.003433171 2.666667e-05
## geometry
## 1 MULTIPOLYGON (((95734.59 12...
## 2 MULTIPOLYGON (((639491.9 22...
## 3 MULTIPOLYGON (((613093.3 85...
## 4 MULTIPOLYGON (((819702.8 20...
## 5 MULTIPOLYGON (((259336 1188...
## 6 MULTIPOLYGON (((565941.6 17...# Map EmissionAve
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties.3435) + tm_polygons(col = "EmissionAve", n=10,
style = "jenks", title = "Emission Ave") +
tm_layout(frame = FALSE, legend.outside = TRUE) 
# Generate 5-mile buffers. Note: 1 mile = 5280 ft (using EPSG 3435)
emission_buffers <- st_buffer(emissions.3435, 5280*5)
glimpse(emission_buffers) #Take a look at the transposed data: check https://dplyr.tidyverse.org/reference/glimpse.html## Rows: 850
## Columns: 14
## $ pointID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ SITE..NAME <chr> "ADM Animal Nutrition", "Prince Agri Products Inc", "…
## $ EIS.Facility.ID <int> 3344211, 3345811, 7314711, 3346911, 3347011, 3362711,…
## $ STATE <chr> "Illinois", "Illinois", "Illinois", "Illinois", "Illi…
## $ State.County <chr> "IL - Adams", "IL - Adams", "IL - Adams", "IL - Adams…
## $ Pollutant <chr> "Lead Compounds", "Lead Compounds", "Lead Compounds",…
## $ Pollutant.Type <chr> "CAP/HAP", "CAP/HAP", "CAP/HAP", "CAP/HAP", "CAP/HAP"…
## $ Emissions..Tons. <dbl> 0.00001, 0.00006, 0.00256, 0.03723, 0.00032, 0.00002,…
## $ Facility.Type <chr> "Unspecified", "Unspecified", "Unspecified", "Airport…
## $ Street.Address <chr> "1000 N 30th St", "4618 Gardner Expy", "401 N Prince …
## $ NAICS <chr> "Other Animal Food Manufacturing", "Other Animal Food…
## $ EPA.Region <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ FIPS <int> 17001, 17001, 17001, 17001, 17003, 17003, 17003, 1700…
## $ geometry <POLYGON [US_survey_foot]> POLYGON ((160215.6 1208455,..., …# Map the buffers
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties) + tm_polygons(alpha = 0.5) +
tm_shape(emission_buffers) + tm_borders(col = "tomato1") +
tm_shape(emissions.3435) + tm_dots(col = "black")
########################
# CALCULATE BUFFER COUNT BY COUNTY
########################
# First we identify how many times buffers overlap counties.
# Then we use lengths() to find out how many items are present in each county.
counties.3435$bufferCt <- lengths(st_intersects(counties.3435, emission_buffers))
head(counties.3435)## Simple feature collection with 6 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95633.61 ymin: 111619.1 xmax: 883631.1 ymax: 2123580
## Projected CRS: NAD83 / Illinois East (ftUS)
## GEOID STATEFP COUNTYFP COUNTYNS AFFGEOID NAME LSAD ALAND
## 1 17001 17 001 00424202 0500000US17001 Adams 06 2214804824
## 2 17003 17 003 00424203 0500000US17003 Alexander 06 609996948
## 3 17005 17 005 00424204 0500000US17005 Bond 06 985073312
## 4 17007 17 007 00424205 0500000US17007 Boone 06 727171389
## 5 17009 17 009 00424206 0500000US17009 Brown 06 791704469
## 6 17011 17 011 00424207 0500000US17011 Bureau 06 2250916543
## AWATER State County LengthLife LLRank TotEmisFac
## 1 41767689 Illinois Adams -0.10 46 4
## 2 44237171 Illinois Alexander 1.50 102 3
## 3 6462629 Illinois Bond -0.13 40 4
## 4 3361806 Illinois Boone -0.40 18 5
## 5 4139310 Illinois Brown -0.80 4 2
## 6 11606402 Illinois Bureau -0.16 38 3
## county_area density EmissionAve
## 1 24325392002 [US_survey_foot^2] 0.004584247 9.965000e-03
## 2 7044304829 [US_survey_foot^2] 0.011872740 6.516667e-03
## 3 10673939103 [US_survey_foot^2] 0.010447277 1.269750e-02
## 4 7864025133 [US_survey_foot^2] 0.017725274 2.994800e-02
## 5 8576601423 [US_survey_foot^2] 0.006501037 3.390000e-03
## 6 24360918399 [US_survey_foot^2] 0.003433171 2.666667e-05
## geometry bufferCt
## 1 MULTIPOLYGON (((95734.59 12... 4
## 2 MULTIPOLYGON (((639491.9 22... 3
## 3 MULTIPOLYGON (((613093.3 85... 7
## 4 MULTIPOLYGON (((819702.8 20... 19
## 5 MULTIPOLYGON (((259336 1188... 3
## 6 MULTIPOLYGON (((565941.6 17... 19tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties.3435) + tm_polygons("bufferCt", n=10,
style = "jenks", title = "CountBuffer") +
tm_layout(frame = FALSE, legend.outside = TRUE)
# Q3: Compare the buffer count by county vs the point count by county.
What difference do you notice? Answer: The buffer count by country is
always much higher than the point count by county. Although the an
emission point is only counted once per exact location, the 5-mile
buffer around each point can interserct multiple counties. Hence,
leading to a much higher count of affected counties when using buffer
zones compared to just using the point locations.
Q4: Further quantify the effects of emission diffusion/dispersion
Integrate what you have learned thus far to answer the following two questions. There are several ways to solve each of the questions.
Q4.1. Calculate the number of counties affected by each of the emission buffer zones.
Write your code chunk below and map the result as a GRADUATED SYMBOLOGY MAP. Be creative in the mapping design.
emission_buffers$affectedCounties <- lengths(st_intersects(emission_buffers, counties.3435))
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(emission_buffers) +
tm_polygons(col = "affectedCounties",
palette = "YlOrRd", style = "quantile", title = "Affected Counties",
legend.is.portrait = TRUE) +
tm_layout(frame = FALSE, legend.outside = TRUE)
Q4.2. Calculate the percentage of area in each county affected by emission diffusion.
Write your code chunk below and map the result as a choropleth MAP.
intersection <- st_intersection(counties.3435, emission_buffers)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriespctAffected <- aggregate(as.numeric(st_area(intersection)), by = list(intersection$GEOID), FUN = sum)
names(pctAffected) <- c("GEOID", "affectedArea")
counties.3435 <- merge(counties.3435, pctAffected, by = "GEOID", all.x = TRUE)
counties.3435$pctAffected <- counties.3435$affectedArea / as.numeric(counties.3435$county_area) * 100
tm_shape(il) + tm_borders(lwd = 2) +
tm_shape(counties.3435) + tm_polygons(col = "pctAffected", n = 10,
style = "jenks", title = "% Area Affected",
legend.is.portrait = TRUE) +
tm_layout(frame = FALSE, legend.outside = TRUE)